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ABSTRACT 


Provided in this study are predicted in-depth temperature and pyrolysis gas pressure distributions for carbon 
phenolic materials that are externally heated with a laser source. Governing equations, numerical techniques and 
comparisons to measured temperature data are also presented. Surface thermochemical conditions were determined 
using the Aerotherm Chemical Equilibrium (ACE) program. Surface heating simulation used facility calibrated 
radiative and convective flux levels. Temperatures and pyrolysis gas pressures are predicted using an upgraded form of 
the SINDA/CMA program that was developed by NASA during the Solid Propulsion Integrity Program (SPIP). Multi- 
specie mass balance, tracking of condensable vapors, high heat rate kinetics, real gas compressibility and reduced 
mixture viscosity’s have been added to the algorithm. In general, surface and in-depth temperature comparisons are 
very good. Specie partial pressures calculations show that a saturated water- vapor mixture is the main contributor to 
peak in-depth total pressure. Further, for most of the cases studied, the water-vapor mixture is driven near the critical 
point and is believed to significantly increase the local heat capacity of the composite material. This phenomenon if 
not accounted for in analysis models may lead to an over prediction in temperature response in charring regions of the 
material. 


NOMENCLATURE 


A -area 

B’ -dimensionless mass loss rate 
C H -Stanton Number, Heat Transfer 
C M -Stanton Number, Mass Transfer 
C -specific heat 
E -activation energy 
F - 1 st generic coefficient 

G -2 nd generic coefficient 

h -enthalpy 
J -mass source/sink rate 
k -thermal conductivity 
K -permeability 
m -mass 
hf -mass flow rate 
mf -mass fraction 
n -number of reactions 
M -Molecular weight 
P -total pressure 
Q -heat transfer rate 
R -gas constant 
¥ -recession rate 
S -source term 
t -time 
T -temperature 
u -velocity 
V -volume 
w -weight fraction 

x -spatial coordinate 

z -compressibility factor 


Z -diffusional driving potential 
a -surface total absorptivity 
P -pre-exponential factor 
e -surface total emissivity 
T -resin volume fraction 
p -dynamic viscosity 
(j) -porosity 
p -density 

y -coefficient for Forchiemer extension 

Subscripts: 
c -carbon 
cn -condensation 
e -edge 
f -final 
g -gas 

i.j -free indices 
o -original 
p -constant pressure 
r -recovery 
rad -radiation 
s -solid material 
sc -solid conduction 
t -total 
v -virgin 
vp -vaporization 
w -wall 



INTRODUCTION 


The Space Shuttle Reusable Solid Rocket Motors (RSRM) have now provided the main propulsion source for 
over 95 missions. During this time a very extensive database of motor operational performance has been amassed 
which includes parameters such as nozzle insulation erosion rates. These data are understood statistically to the extent 
that variations on the order of tenths of an inch are indicators that a change has occurred in either materials and/or 
processes used in their refurbishment. The left hand nozzle of the RSRM-56 flight set displayed anomalous erosion 
(pocketing) aft of the throat (Fig. 1) affecting the full circumference of the motor and measuring as much as 0.5” deeper 
than expected mean values. Based on statistics, the event was approximately a 6-cr occurrence and thus could not be 
discarded without further understanding. 
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Fig. 1. RSRM Nozzle Showing Pocketing Region 

The degraded material performance was believed to be attributable to the “pocketing” phenomenon that is 
distinctly different from typically occurring thermochemical erosion. At this location in the nozzle throat ring material 
ply angles are 45° to motor centerline and about 70° to the conducted isotherms. It is known that in-plane (with ply) 
fibers oriented orthogonal to the isotherms are more likely to pocket. It was therefore suspected that for the RSRM-56 
nozzle, process variation had produced fiber orientations approaching 90° to the flame surface and was likely the 
primary cause of the increased erosion. Additionally, other factors related to materials and/or process variation were 
considered potential contributors thus it was decided to initiate a comprehensive test program aimed at gaining a better 
understanding of material thermostructural behavior. 

The resources of the Laser Hardened Material Evaluation Laboratory (LHMEL) facility were utilized to examine 
pocketing activity as a function of fiber orientation and other material variations such as resin content, moisture content 
and ply distortions. LHMEL has the major advantages of a relatively large spatially flat surface heating distribution of 
precise magnitude, rapid turn-around test time and direct measurement of surface temperature. Disadvantages of the 
LHMEL are total pressures, thermochemistry and surface recession does not compare well with the actual RSRM. 
Average recession rates are about one-forth of that experienced in the RSRM nozzle at the location of interest. There is 
some debate and conflicting data ^ that seems to suggest that the effect of active surface thermochemistry may be 
important in terms of suppression of pocketing. Notwithstanding these data, the decision was made to test at LHMEL 
based on the belief that pocketing is an “in-depth” phenomena and not strongly dependent on surface recession. 

The following provides a description of modifications incorporated into the S1NDA/CMA computer code which 
was developed by the author ^ during the Solid Propulsion Integrity Program (SPIP). Upgrades include multi-specie 
mass balance, real gas equation of state using generalized compressibility data, reduced mixture viscosity, resin 
weight fraction Arrhenius formulation, high rate TGA coefficients and a condensation/ vaporization simulation for 
vapors in the pyrolysis gas mixture. Basic formulations of the energy and momentum equations remain essentially 
unchanged but will be covered for completeness. 



Governing Equations 


In-depth temperature and pyrolysis gas pressure calculations are based on simultaneous solution of 1-D 
conservation equations for mass, momentum and energy along with a real gas equation of state and kinetic rate 
equation. The following basic assumptions were incorporated into the math model: 

1) Local thermal equilibrium exists between pyrolysis gas and solid thus one energy equation can 
describe thermal response of both. 

2) Pyrolysis gas motion is governed by the Darcy-Forcheimer equation. Permeability and porosity data 
was correlated as a function of degree of char. Data was assembled from Clayton^ and Stokes^. 

3) Temperature and pressure gradients are 1-dimensional thus material anisotropy can be simulated by 
use of effective properties. Rule of mixtures was used for determination of properties in the charring 
region. 

4) Transport of condensable species through the pore network occurs in the vapor phase. Liquid 
occupying pore volumes was assumed to be stationary' and in equilibrium with its respective vapor in 
the mixture. Condensation and vaporization rates are governed by the amount of a specie that can be 
thermodynamically accommodated in the mixture relative to saturation over a given time step. 

5) Condensed phase species residing together in a pore volume are assumed not to interact with each 
other chemically or physically. Gas phase permeability remains unchanged due to the presence of 
liquid in the pore volumes. 

6) Mixture specie concentrations, in the pore free volumes, were determined by “origin” generation, 
condensation/vaporization rates and upstream advection. Equilibrium and/or kinetic rate reactions 
within the gas and reactions with the char layer are currently not modeled. 


Surface Energy Balance 

Surface heating conditions are determined by consideration of combined convective, radiative and 
thermochemical loading. The LHMEL is unique, and different from the RSRM in that surface response is driven by 
the incident radiation emitted from a C0 2 laser and is convectivly cooled by air flow. Oxidation of carbon in the char 
layer is present and averages about one mil/sec depending on incident heat rate and location on the sample. 
Components of the surface heat flux are depicted in Fig. 2. 
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Fig. 2. Thermochemically Eroding Surface Boundary Conditions 

During testing, surface radiometer data was collected and used for model calibration of surface optical properties. 
Backfitting model response to measured data across the range of incident flux levels, surface absorptivity was found to 
be independent of temperature while re-radiated energy levels were controlled by a temperature dependent emissivity. 
The backlit suggested that a ~ 0.97 and e was ~ 0.85 @ 3000°F and increased linearly to ~ 0.96 @ 5000°F. Forced 
convective cooling of the sample was imposed by a 0.5 Mach air flow directed parallel to the heated surface. Facility 
airflow calibration data was used for determination of convective heat and mass transfer coefficients and average values 
of these quantities input into the SINDA/CMA model. Surface oxidation rate was correlated in familiar b-prime table 
format and estimated by an ACE ^ solution for standard air environment. Surface recession rates average about one 



mil/sec and vary substantially along the heated surface with the leading edge having the greatest erosion. Arithmetic 
summation of respective fluxes identified in Fig. 2. gives the following expression for the surface energy balance: 
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Terms left to right are identified as free stream convection, surface thermochemical, backside advective fluxes, 
radiation incident/emitted and conduction into the material. Numerical solution for the surface energy balance^ is 
implicit with respect to temperature calculations but explicit in time relative to mass flow calculated quantities. Net 
flux values are loaded into the source term of the surface node during iterative convergence of the global temperature 
and pressure calculations. This numerical method is different from the standard CMA approach but the two methods 
compare well with a difference of less than 1/2 % in computed surface temperature^. 


In-Depth Thermal Solution 

Invoking the assumption of gas-solid equilibrium, leads to the standard CMA^ formulation for energy 
conservation given by Eq. 2. For this study an additional term that accounts for pyrolysis gas capacitance has been 
added to the equation per the general formulation provided by Keyhani ^ l Inclusion of gas capacitance has heretofore 
been considered unnecessary due to order of magnitude considerations however based on findings presented in this 
study, it is believed this term can become significant in charring regions of the material. Terms left to right are energy 
storage, conduction, decomposition, grid movement, pyrolysis gas flow and latent phase change rates. The last two 
terms were added to account for the phase energy of saturated water and phenol compounds. 
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A finite element scheme was used for discretization of the energy equation. The computational grid consist of 
one dimensional first order elements with applied front/back face boundary conditions, Fig. 3. 
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Fig. 3. One Dimensional Finite Element Grid 

Surface recession is accommodated by movement of the grid relative to a fixed coordinate and applying a 
correction term to account for the induced adverted energy into the elementfs). As with the baseline CMA technique, 
the last element (furthest from the heated surface) shrinks to accommodate surface recession. If eroded total depth 
exceeds elemental thickness it is dropped from the active network and the adjacent element now absorbs the recession 
and so on. Determination of the elemental “stiffness” matrix is based on trapezoidal rule numerical quadrature which 
evaluates material heterogeneous and temperature dependencies explicitly in time. Temperature and pressure elemental 
integration points are coincident with a density field “nodlet” grid. Use of a nodlet grid for the density calculations is 
similar to the CMA technique and is generally necessitated by the exponential behavior of the Arrhenuis equation. 




In-Depth Pyrolysis Gas Pressure Solution 


Pyrolysis gas pressure distributions are calculated using a Darcy-Forcheimer form of the momentum equation 
substituted into the conservation of mass equation. Real gas effects for the mixture are simulated by application of a 
generalized compressibility factor to the ideal gas law. Pseudocritical temperatures and pressures^ are calculated for 
the mixture based on mole fractions and the individual specie data. These “reduced” properties are used as 
independent variables for table lookup to determine the z factors. Mixture viscosity calculations^ incorporate the 
effects of pressure, temperature and molecular polarity and are functions of mole fraction data and the pseudocritical 
reduced properties. Expression of gas properties in terms of mixture equivalents, i.e., gas constant, viscosity’s, etc., 
permits use of Darcys equation for computation of total gas pressures. Specie partial pressures are simple functions of 
the mole fraction data which result from the multi-component mass balance. Terms left to right in Eq. (3) are 
functionally identified as gaseous mixture storage, total diffusive mass flux, total rate of pyrolysis gas generation, 
coordinate set movement correction, total multi-specie vaporization and condensation. 
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Eq. (3), the “pressure” equation, is analogous in formulation to the energy equation and is thus numerically solvable by 
the same finite element techniques used for temperature calculation. Assembly of the pressure elemental stiflhess 
matrix relies on explicit quadrature of spatially dependent properties using the same procedure derived for the thermal 
calculations. Details of discretization of the diffusive term, treatment of source/sink terms and numerical degree of 
coupling are developed and discussed in Ref. [2]. 


In-Depth Kinetic Decomposition 

A resin based weight fraction formulation of the Arrhenius equation was employed for the in-depth material 
decomposition simulation. Kinetics coefficients were developed by Clayton^ 10 ' and computed from high rate TGA data 
derived by Southern Research Institute (SoRI). The weight loss curve fit considered three reactions and was based on 
the 3000°C/minute data. The Arrhenius relationship used in the SINDA/CMA code has the following form: 
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A thermal decomposition “nodlet” grid was implemented for refined description of variation in material density 
along the 1-D coordinate. This grid is fixed and contains time based composite material density resulting from 
integration of Eq. (4). A simple Euler scheme was applied in which integration time steps are set equal to the transient 
solution time step for the temperature and pressure calculations. Updating of temperature in the calculation occurs 
explicitly and the resultant pyrolysis gas flow rates are used explicitly in temperature and pressure calculations (loosely 
coupled). Quadrature involved in evaluation of stiffness matrices assumes a piecewise linear distribution of density 
described by integration of Eq. (4) at the fixed nodlet sites. Conversion of resin weight fraction data to composite 
density was based on the following relationship 
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Multi-Component Mass Balance 

A multi-component mass balance allows for tracking of individual molecular species evolved during the in-depth 
thermal decomposition process. The procedure utilizes a control volume aligned with elements in the pressure grid, 
Fig. 4. Total rate of decomposition and thus pyrolysis gas generation is determined by Eq. (4) and is assumed strictly a 
function of temperature and local char state. Mole fraction data that describes molecular species evolved as a function 



of degree of char taken from Clayton^ was used to determine the “origin” generation rate. A total of eight molecular 
species were considered; water, carbon dioxide, carbon monoxide, methane, hydrogen, phenol, cresol and xylenol. Of 
these eight, water and phenol have critical temperatures high enough and occur in sufficient concentrations that 
condensation and vaporization has to be considered if accurate total pressure magnitudes are to be calculated. 



At the boundaries of the control volume are the advected fluxes of the individual species. Depending on computed 
direction of flow upstream contributions contribute to a weighted average type calculation of pore volume specie 
concentration. Equivalent molecular weights, gas constants, mixture viscosity’s, specific heats are all functions of the 
mixture mole fraction calculations. Tracking partial pressures of individual species allows for simulation of 
condensation and vaporization. The computational procedure involves comparing specie partial pressure to its 
saturation pressure for the local temperature. Fig. 5. Below critical temperatures for the given specie, if its partial 
pressure tries to exceed the saturation pressure, an instantaneous rate of condensation is calculated that will keep the 
specie partial pressures equal to its saturation pressure (Tj to T 3 ). Time integration of the rate of condensation gives 
the total amount of liquid that has accumulated in the open pore volumes. This liquid is available for vaporization 
when conditions are such that the mixture can thermodynamically accommodate its presence. Vaporization rates are 
computed based on the premise that the gas mixture remains saturated until all the liquid in the pore volume is 
removed(T 5 to T 5 ). Above critical temperatures gas mixture PVT behavior is described by the ideal gas law using 
generalized compressibility factors. The technique employed is similar to that used by Clayton^. 
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Fig. 5. Condensation and Vaporization Simulation 
NUMERICAL SOLUTION 

As previously discussed, discretization of energy and mass conservation equations was based on a finite element 
formulation employing a nodlet grid for integration of the Arrhenius equation. Time integration of the non-steady 
behavior of the diffusion equations, i.e., pressure and temperature, was performed by a Crank-Nicholson procedure. 
Eqs. (2) and (3), are cast into the following generic form: 
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where; O is temperature or total pressure 


Eq. (5) is solved by iteration for the dependent variable AO ; . A successive point (Gauss-Seidel) scheme is 
applied which uses dependent and source term damping. Coupling is folly implicit between pressure, temperature and 
surface energy balance meaning that all quantities are converged together along with their respective coefficients and 
source terms at every time step. The Arrhenius equation is not iterated with temperature and pressures. A simple 
explicit updating is performed using “old” time step data values. The global method of solution is described in detail 
by Clayton^ but will be summarized here to include the phase change logic. The overall numerical procedure goes as 
follows: 

1) Initialize temperature, pressure, density and nodal coordinates. Compute coefficients in Equation (5). 

2) Increment boundary information and solve for temperatures and pressures by iteration. 

3) Using converged data in Variables 2, interpolate temperatures onto density grid and integrate Arrhenuis 
equation across the time step. Store decomposition data into an array versus position. 

4) Recalculate coefficients in Eqn. (5) based on new properties data, i.e., conductivity's, permeability's, mass- 
energy source and sink rates, coordinate system location, etc... 

5) Perform nodal mass balance as function of converged flow conditions for current time step. Compute mass 
and mole fractions, partial pressures, mixture equivalent properties. 

6) Compare partial pressures with saturation pressure @ temperature for condensable species. If partial pressure 
is greater, using real gas law compute amount of mass removal necessary to make the two equal. Accumulate 
this mass as liquid in the control volume. If partial pressure is less than saturation, using real gas law compute 
amount of mass necessary to saturate mixture and vaporize accumulated liquid (if there is any). Adjust mass 
source terms in pressure network to reflect local rates of condensation/vaporization. 

7) Perform grid movement logistics, if current time is less than end time return to step #2 

Steps #3-#7 are performed in Variables 2 of the SINDA/CMA model thus all procedures described in these steps are 
explicitly coupled in time to the pressure and temperature calculations in step #2. 


RESULTS 

Spatial distributions for in-depth temperature and pyrolysis gas pressure for various incident radiant heating rates 
are presented in Figs. 6- 14. For these cases ply angles are fixed at 90° and time slices at 3, 10 and 20 seconds are 
provided. For clarity, partial pressures of only major contributors such as water vapor, carbon dioxide and monoxide 
are presented. For the 300 Watt case, surface temperatures range from 2900°-3700°F and increase monatonically 
during the test. Peak total pressures range from 140-180 atmospheres with the maximum occurring at the 10 second 
time slice. Water vapor is the dominant pyrolysis gas specie in the charring regions of the material with mixture mole 
fractions approaching 99%. Gas flows are driven in-depth and to the heated surface depending on proximity relative to 
the peak pressure location. In the cooler material, carbon dioxide becomes the dominant specie. As the mixture is 
driven in-depth, the water vapor is condensed out leaving only species with critical temperatures low enough to exist as 
a gas at the given total pressure and temperature. For the 500 Watt cases. Figs. 9-11, surface temperatures range from 
3600°-4300°F, increasing during the test. Peak total pressures range from 130-190 atmospheres with the maximum 
occurring at 3 seconds (earlier in test compared to 300 Watts). Specie distributions follow the same general trends. 
Looking at the 1000 Watt cases. Figs. 12-14, surface temperatures now vary from 4600°-5 100°F and increase during 
the test. Peak total pressures range from 200-220 atmospheres with the maximum occurring at 10 seconds. Clearly 
observable trends in Figs. 6-14 are overall increase in material temperatures and total pressures with increasing surface 
heat flux. 



temperature, °F 



Fig. 6. Temperature and Pressure Distributions @ 3 Seconds 
300 Watt Case, 90° Ply 
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Fig. 11. Temperature and Pressure Distributions @ 20 Seconds 
500 Watt Case, 90° Ply 
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Fig. 12. Temperature and Pressure Distributions @ 3 Seconds 
1000 Watt Case, 90° Ply 
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Fig. 13. Temperature and Pressure Distributions @ 10 Seconds 
1000 Watt Case, 90° Ply 
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Fig. 14. Temperature and Pressure Distributions @ 20 Seconds 
1000 Watt Case. 90° Ply 


Provided in Figs. 15 and 16 are comparisons of predicted material thermal response to the laboratory measured 
data. Surface temperature data was obtained by radiometry while in-depth measurements were gathered by secondarily 
bonded thermocouples. Thermocouple depths were 0. 125. 0.250 and 0.375 inches; a constant spacing of one-eighth 
inch. Adjustment of surface optical properties, per the discussed procedure, allowed for very good correlation in 
predicted response to the measured data. A good temperature match at the surface is a necessary starting point for 
understanding comparisons made at the in-depth locations. In general, the in-depth calculated quantities compare well 
with measured data and some discernible trends were evident. For the 300 Watt case at the 0.25” location, the 
measurement strays from the prediction at —10 seconds. Comparisons at locations on both sides of this thermocouple 
are very good thus it is believed that the measurement may be inaccurate to some extent. High surface heating rate 
(1000 Watt case) comparisons are provided in Fig. 16. Unfortunately, lower capability thermocouples were used for 
this test and the in-depth measurements fail at — 2100°F. A general tendency for this test, and others in a series of 
-1100 tests, is that model predictions lead measurements in the 700-1500°F range. Laboratory measured specific 
heats, for virgin and char material, were used for results presented in Figs. 15 and 16 to demonstrate this over 
prediction trend. Addition of saturated water pressure and temperature dependent specific heats into the energy 
equation would produce a better temperature comparison at the in-depth locations. 








CONCLUSIONS 


Based on findings presented in this study the following conclusions are made: 

1) Surface thermal simulation was best backfit by assuming a constant absorptivity and temperature 
dependent emissivity. At surface temperatures approaching 6000°R the two are equal at ~ 0.97 while 
at the lower temperatures, emissivity values were estimated to be ~ 0.85. It is recognized these 
backfit values are sensitive to the assumed radiometer values used during testing. 

2) In-depth thermal response is not strongly dependent on detail calculation of the pyrolysis gas flow 
field. Somewhat satisfactory results have been obtained for years assuming gas flow is always 
directed to the heated surface and vapor condensation not a factor. The reason for the “weak” 
coupling is that in-depth thermal response is driven primarily by conduction into the material. 
Pyrolysis gas flow contributions to the overall energy balance are second order effects. 

3) The trend of increasing total pressures with increasing surface heat rate is attributable to material 
“kinetic shift” meaning basically that at the higher heat rates, the material has a tendency to be less 
charred at higher temperatures. Trapped volatile's and initially evolved gases are dealing with higher 
temperatures and logarithmically smaller shifts in permeability thus pressure build up is greater. 

4) Not accounting for pyrolysis gas reactions with carbon in the char layer seems to be a reasonable 
approximation at temperatures < 2000-2500 :, F. This premise is supported by findings presented by 
April^ 111 were specie concentration data was obtained for gas flow through char layers at various 
temperatures. Peak magnitudes of pyrolysis gas pressure build up, see Figs. 6-14, take place in 
partially decomposed material where local temperatures are in the 700-1 100°F range. Water-carbon 
reactions within the char layer could potentially increase local permeability and thus affect pressure 
magnitude and distribution obtained from the global solution. The exact extent of influence is 
unknown at this time and suggest that permeability may be correlated versus actual material density 
rather than the degree of char parameter. This method of correlation could potentially capture the 
effect of residual char density changes due to heat rate dependence and/or enhanced pyrolysis gas 
reactions with carbon. 

5) For a given heat flux, calculated gas pressures for ply angles less than 90° are greater than pressures 
calculated for the 90° case. This is a result of the across-ply permeability component coming into 
play in the effective 1-D property calculations, i.e.. across-ply « in-plane permeability’s at 
temperatures less than ~ 750°F. Gas generation rate is essentially unchanged while flow resistance 
has increased thus in-depth pressure build-up is greater. This trend is based on the premise that 
permeability is a function of degree of char only which is how the data was correlated in the thermal 
model. Its is known that permeability can be a function of compressive load which has the 
implications that the overall solution will necessarily have to couple thermal and structural response. 

6) Formulation of the energy equation includes the local heat capacity of pyrolysis gas as contributing to 
the storage of energy in the material. The advective terms have always been included in CMA type 
codes but storage terms neglected on the premise of being second order. Results provided by the 
multi-specie calculations indicate that a liquid water-vapor mixture can exist during the 
decomposition process and that the mixture can be driven near critical conditions. In theory, a 
substance at the critical point has an infinite heat capacitance^ 2 ^ and the asymptotes, near the 
singularity, are finite and are thermodynamically obtainable to a fixed extent. Historically, there has 
been a tendency to over predict in-depth temperature response using laboratory measured thermal 
properties. Many theories have been proposed to explain the differences which include kinetics, 
dynamic conductivity’s, instrumentation, but it is believed by findings presented herein that part of 
the in -accuracy may be a result of not considering the thermodynamic state of water and implications 
of its pressure and temperature history. 
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